function X = h2v(N)
% If A is a symmetric n-by-n matrix, then let vech(A) denote the n*(n+1)/2
% vector obtained from vec(A) by eliminating all of the supra-diagonal
% elements of A (as in Magnus and Neudecker). The duplication matrix D is: 
% 
%    D*vech(A) = vec(A) 
% 
% Written by Alan Sutherland and modified by Michael Johnston
nr = N^2 ;
nc = N*(N+1)/2 ;

X=zeros(nr,nc);
nn=2*nc-nr;
ir=1;
ic=1;
for ii=1:nr;
    if ir<=ic
        jj=(ic*(ic-1)/2)+ir;
    else
        icc=ir;
        irr=ic;
        jj=(icc*(icc-1)/2)+irr;
    end
    X(ii,jj)=1;
    ir=ir+1;
    if ir>nn
        ic=ic+1;
        ir=1;
    end
end
end